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I. INTRODUCTION 

The mysterious pseudogap (PG) phase is one of 
the most fascinating aspects of the underdoped high- 
temperature curpate superconductors (HTCS). By a va- 
riety of probes the pseudogap (the suppression of the 
low-energy single-particle spectral weight) has been ob- 
served to persist from above the superconducting (SC) 
critical temperature T c to T* in the underdoped regime. 
The direct evidences of this spectral gap come from the 
angle-resolved photoemission spectroscopy ( ARPES) EES. 
Ding et al\^ studied the underdoped Bi2Sr2CaCu20s+5 
using ARPES and found that a pseudogap with d-wave 
symmetry begins to open up for T < T* and develops 
smoothly into the d-wave SC gap below T c . One pecu- 
liar property of the pseudogap phase revealed by further 
experimental investigation is the truncated Fermi surface 
termed as Fermi arcsp exhibiting distinct difference from 
the point-like (four gap nodes) Fermi surface for T well 
below T c as expected for a pure d-wave superconductor 
and the closed Fermi surface for T above T* . 

There are two basic scenarios of the PG phase. The 
first one attributes the opening of the pseudogap to the 
presence of an e xoti c order competing with the SC phase, 
such as the spirl^ and/or charge 7 density waves and so 
on. The second scenario associates the PG phase with 
the phase-incoherent pairing and therefore the pseudo- 
gap is interpreted as a precursor of the SC order. In this 
preformed-pair scenario there are two energy scales: one 
is the BCS energy gap A which is closely related to the 
binding energy of the electron pair, and the other is the 
phase-stiffness energy scale Tg which protects the phase 
coherence. For the conventional superconductors Tg is 
larger than A so that the SC state is destroyed by pair 
breaking. However for the underdoped HTCS, because of 
the low carrier density and the short correlation length, 
A is larger than Tg^ and therefore the phase coherence 
is destroyed while the energy gap survives as tempera- 



ture increases across T c . In this context T c is determined 
by Tg and the pseudogap is caused by the pair fluctua- 
tion^sHlHl Franz and Millis^ showed that random super- 
current induced by thermal phase fluctuations can cause 
the shift of electronic spectral weight in both momentum 
and energy. Berg and AltmarPSl further attributed the 
emergence of the Fermi arc to the pile up of the low- 
energy spectral weight along the underlying Fermi sur- 
face due to the Doppler-shift effect of the fluctuating su- 
percurrent. This picture of phase fluctuations is concise 
and instructive, yet the analytical results relied on the 
semiclassical approximation^ where only far-field effect 
of the vortex-type excitations is considered, which might 
be uncontrolled as argued in Kef. Hi . Furthermore, the 
probability distribution of the fluctuating supercurrent 
was assum ed ph enomenologically to be Gaussian type. 
Recently w o ' 17 ' * 18 ! attempted to go beyond the semiclassi- 
cal approximation by employing a 2D XY model to sim- 
ulate the vortex-type phase fluctuations and numerically 
taking both the Doppler effect of the whirling supercur- 
rent and the scattering effect of vortices as topological 
singularities into full consideration. However, the XY 
model is still a phenomenological description of the phase 
fluctuations, which includes a temperature-independent 
phase-stiffness constant J. 

In this work, we start from a 2D attractive Hubbard 
model with only nearest-neighbor interactions to investi- 
gate the pseudogap phase and the evolution of Fermi arcs 
in d-wave superconductors. The path-integral formalism 
is employed where pairing fluctuations are inherently em- 
bedded. A local-update Monte Carlo scheme on the basis 
of the Green's function method is presented to speed up 
the random walk in the classical configuration space of 
pairing field. Superfluid density is calculated as the sig- 
nature of the SC phase transition and compared with the 
phase correlation function^. The single-electron spec- 
tral function is calculated using Chebyshev polynomial 
methocP^H^l The temperature dependence of Fermi-arc 
length is discussed. 



2 



The paper is organized as follows: In Section [TT] we de- 
scribe the basic path-integral formalism to treat the 2D 
extended Hubbard model. The local-update algorithm 
based on Green's function theory and the Chebyshev ex- 
pansion approach are presented. In Section |HI] we calcu- 
late the temperature dependencies of various quantities 
relevant to the phase fluctuations and pseudogap phase. 
The conclusion is given in Section |IV| 

II. THE MODEL AND FORMALISM 

The BCS Hamiltonian Hbcs we adopt is given by: 

#BCS = ~t ^ 4a C t+Sa ~ *' Y C L C *+5V 
z,(5, cr i.S' ,er 

i,5 i.cr 



Here <r denotes spin, i is the index of site of the two- 
dimensional L x L square lattice, i + S and i + 8' denote 
the nearest-neighboring (NN) and next-NN sites of i, re- 
spectively t and t' are the NN and next-NN hopping 
integrals, respectively, /j. is the chemical potential. The 
attractive interaction V > between electrons on the 
NN sites favors unconventional superconducting phases. 
To investigate the effect of superconducting fluctuations, 
the quantum partition function is expressed in the path 
integral formalism l 21 l 22 l 



Z= D{<p irT (T),$ i(T (T)}exp(-S) (2) 



where the action S is expressed as 



S((p,cp) 



ch 



tijip l!y {T)(f J<y (T) — V22 ( Pit{ r ) t Pi+Sl( T ) ( Pi+Sl( T ) i Pit( T ) 



(3) 



where ip^ and (f>i a denote Grassmann fields and j3 — 
l/ksT. We then decouple the quartic term in the action 
by introducing an auxiliary Hubbard-Stratonovich field 
Ai ! i + 5(r) in the Cooper channel. For a square lattice 
with TV sites and periodic boundary condition, there are 
totally 2N(N = L 2 ) independent A iji+( 5(r)'s. Hereafter 
we use A to denote the set {A^+^r)}. The partition 
function now becomes 

Z= f DADAe-^ A > A \ (4) 

where 

J DADA = J HYl dA iii+g (r)dA iti+s (T). (5) 

i—1 S—x,y 

In Eq. Q, the grand potential is expressed as 

Q(A, A) = Q,(A, A) + V- 1 Y |A m+5 (t)| 2 , (6) 

i,5 

where Q f denotes the fermionic thermodynamic potential 

%(A,A) = -/3- 1 lnTre~' 3dBdG ( A ). (7) 
Here the BdG Hamiltonian is written by 

i?Bd G (A) = ¥H BdG (A)V (8) 




where W (<I>) denotes the Nambu creation (annihilation) 
operator defined as * f = (c| t , c n , c| t , c 2 j., • • • , cj^, c Ni ). 

&BdG is a 2 N x 2 N Hermitian matrix which will be called 
BdG matrix. Hereafter we use capital letters with a tilde 
( ~ ) to denote 2N x 2N dimensional matrices (e.g. the 
BdG matrix -ffedc) while a hat (") to denote operators in 
second quantization (e.g. the BdG Hamiltonian Hsdc)- 
For the sake of convenience, we will omit the argument 
A and simply use -ffedG an d i?BdG to denote the BdG 
Hamiltonian and matrix for a certain pairing held A. 

In the following, we will ignore the T-dependence of 
A ! ; i i+a(r), i.e. the quantum fluctuation, and concentrate 
on its thermal fluctuations expected to be dominant near 
T c especially in the high temperature pseudogap region. 
With this approximation, i/edG actually describes the 
electrons moving in a static but spatially fluctuating pair- 
ing field. Moreover Eq. Q becomes a classical partition 
function expressed as an integration over the classical 
phase space formed by {Ai^ + s = \ Ai : i + s\e l '^ i }, whose di- 
mension is 47V for a TV-site square lattice. Such multidi- 
mensional integration can be performed by the standard 
Monte Carlo method. To achieve this goal, one need to 
obtain the probability distribution P(A) cx e~^ A ' A ) 
for a configuration A, or its change characterized by the 
ratio 



-/3[n(A',A')-n(A,A)] 

P(A) 
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for a possible change of configuration A — > A'. The 
acceptance probability for such a change is given by 



P A (A' «- A) = min 



P(A) 



(10) 



according to the Metropolis algorithm. To obtain es- 
pecially the nontrivial 51/, previous numerical worl^l^l re- 
lated it to the eigen-spectrum of the BdG matrix through 
the relation S1/(A, A) = -/3" 1 J2 n M 1 + e" 13 ^), where 
e n is the eigenvalue of the BdG equations, 



E 
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11? 
u n 

1)0 



(11) 



To solve this eigenvalue problem, one needs to diagonalize 
a 2N x 2 A BdG matrix and the workload is 0(A 3 ) which 
is quite time-consuming for large lattice. 

In this paper, we will propose an alternative local- 
update scheme based on the Green's function method. 
The Gor'kov Green's function is employed, which is de- 
fined as: 



G(iT,jr') 



-T. 



C i r( T )c] T (T') C.4 0)0^(7-') 



(12) 



in terms of 2 x 2 Nambu matrix notation. Here (• ■ • ) = 
Tr[- • •e - ' 9 ^ BdG ]/Tr[e - ^ BdG ], and its Fourier transform 
with respect to the imaginary time is 



G(i,j;iuJk) 



(13) 



where u>k — 2-KT(k + 1) the Matsubara frequencies. One 
can easily find that the matrix form of the Gor'kov 
Green's function is actually the resolvent of the BdG ma- 
trix H B dG, i-e. 



G(iw fe ) = (iu) k I- HBdcT 



(14) 



whe re I the unity matrix. Combining this relation with 
Eq. ( 11 1, we obtain the spectral representation of G(iwfe), 



E 



lUJ k 



(15) 



which is a 2 x 2 matrix. Before starting our simulation, we 
only need to diagonalize the BdG matrix with certain ini- 
tial (often random) configuration once for all. Then the 
eigen-energies e„ and eigen-functions u n and v n are used 
to calculate the Green's function according to Eq. (15 1. 



As we will show in Section II A we can always update 
the Gor'kov Green's function without having to diago- 
nalize the BdG matrix any more as long as the change of 
configuration is proposed locally. 



A. Update of the Green's function and calculation 
of the acceptance probability 

We assume a local change of the configuration, located 
at the 1 and 2 sites without loss of generality, from Ai^ to 
A[ 2 = A12 +Xi,2- Then the BdG Hamiltonian becomes 



H' 



Hp 



BdG - -"BdG I" Hi 

Hi = Xi,2(4t c 24. + 4f c i4.) + h - c - 



(16) 
(17) 



where Hi denotes the corresponding change of the BdG 
Hamiltonian. According to Eq. ( 14 ) , we have 



&{iu n ) = {tuj - H BdG - Hi)- 1 

= G{iu n )[l - HiG{iuj n )]-\ (18) 

where Hi denotes the matrix form of Hi in the Nambu 
representation as Eq.([8|), 



Hi = 



Xax4 












(19) 



2Nx2N 



where only its upper-left-corner 4x4 block is has non- 
zero elements, and is 







^4x4 — 



Xl,2 \ 





xl,2 
Xi,2 o 








(20) 



The inverse operation on the right hand of Eq. ( 18 ) can 
be performed as follows, 



(I-HiG)- 1 = 



I - 





(I-XA)- 1 (I - XA)~ 1 XB 
/ 



(21) 
(22) 
(23) 



In the above derivation, we use block matrices A, B, C, D 
to denote G, whose dimension is 4 x 4, 4 x (2 A — 4), 
(2A-4) x4 and (2A-4) x (2A-4), respectively. As most 
non-zero elements of the above matrix is concentrated 
on t he fi rst four rows, one can update G' according to 
Eq. (18) with 0(N 2 ) computing operations. 



Next, we will shown how to obtain the change of the 
thermodynamic potential which determines the accep- 
tance of the proposed local update. According to text- 
book^!, one nas 



(24) 



n' f -n f = / d\(Hi) 
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where (•••)* = Tr[- ■ ■ e^^/Trfe-^^)] with H(X) = 
-ffsdG + XH\. The integrated function (Hi)\ are also 
related with the Gor'kov Green's function^ 

(#i) A = 2Re{xi 2 [G A (lr, 2r+) 21 + G A (2r, lr+) 21 ]} 
- 2Re{ X i 2 T^[G A (l, 2; iu k ) 21 + G A (2, 1; lWfc ) 21 ]} 

iuik 

(25) 



Using Eq. (18) and (23), we have 

G A (l,2;iw fc )2i = [G(iuj k ) ixi (I - XX A)- 1 ' 
G A (2,l;iw fc )2i = [G(m.)4x4(/-A1A)- 1 ; 



2..", 



(26) 
(27) 



where |k f) = iV -1 / 2 ^ - e lk l cL|0). Generally, one can 
first solve the BdG equation ( fll| ), then employ the fol- 
lowing equation, 

A A (k, w) = £ 5 (" - ^)«^ MH) (33) 



to calculate the spectral function for one configuration. 
However, since the computational effort of full diagonal- 
ization of the BdG matrix is Q (N 3 )j we will apply the 
Chebyshev polynomial approac h 1 18 ! 19 !, which is O(MN) 
with M -C N 2 , to cut the computational cost. 
We perform a Chebyshev polynomial expansion 



Therefore, from Eqs. (25), (26) and (27), the integration 



over A in Eq. ( 24 ) is readily performed 



S(x - y) 



7T\/l 



I' 



+ 2j2T m {x)T rn {y) 



, (34) 



/ (I-XXA)- 1 dX = -(XA)- 1 ]ji(I-XA) (28) Eq (|32j). we have 
Jo 



to handle the Dirac S function. After substituting it into 



-0 



ln(l-di) 
dl 








\ 





ln(l-da) 
d-i 














ln(l-d 3 ) 
da 














ln(X-d 4 ) 
<2 4 / 



A A (k,a,) = Mo + 2E ^ 1 ^g^ (^) i (3fi) 



where 



/i m = (kt|T m (i?Bd G /s)|kt), 



(36) 



where to treat the 4x4 matrix as an argument of loga- 
rithm function in Eq. ( 28 ), we diagonalize XA = ODO^ 1 



with O the transformation matrix and D the diagonal 
matrix with eigenvalues di,2,3,4- 



B. Spectral function and Chebyshev polynomial 
approach 

With the help of the local-update scheme described 
above, the classical phase space of A are sampled and 
thermodynamic averages of physical quantities can be 
obtained. As an example, we give the definition of the 
single-electron spectral function, 

A/ s [ DADAA A (k,uj)e-^ A ^ , x 

where A A (k, ui) denotes the single-electron spectral func- 
tion for a certain configuration of A. A A (k,uj) can be 
derived according to 



are Chebyshev moments. Here for numerical calculation, 
the infinite series in Eq. ( 34 ) has to be truncated by M as 
shown in Eq. ([351 and to damp the consequential Gibbs 



A A (k,u;) = --ImVG R ^»ne lk(i - 

7T ' * 



(30) 



where the retarded Green's function is the real-space rep- 
resentation of the resolvent, 

G R (z, j,u) n = {% t |(w + i0+ - HBdcr^j t). (31) 
Combining the above two equations, we have 

A A (k > w) = <kt|$(«-.ffBdG)|kt> J (32) 



oscillations the Lorentz kernel g m is used in Eq. ( 35 1 with 
g m = sinh[A(l — m/AT)]/ sinh(A) where A is a free param- 
eter of the kernel and we choose A = 4 throughout our 
calculation as a compromise between good resolution and 
sufficient damping of the Gibbs oscillations as suggested 
in Ref. [19]. s denotes the scaling factor ensuring the 
spectrum of H^q/s falling into the interval [— 1, 1], i.e. 
the domain of the Chebyshev polynomials. 

Most computational effort is spent in the calculation 
of the Chebyshev moments fi m according to Eq. ( |36| , 
which reduces to sparse matrix-vector multiplications af- 
ter taking advantage of the recursion relation T m (x) — 
2xT m _i(x) — T m _ 2 (x). Considering that the BdG Hamil- 
tonian is sparse, the cost of matrix- vector multiplication 
is an O(N) process and the calculation of M moments re- 
quires only O(MN) computational operations. Further 
relations of the Chebyshev polynomials T 2m = 2T^ — 1 
and T 2m +i = 2T m T m+ i — T\ enable us to obtain two mo- 
ments per matrix-vector multiplication. Therefore, cal- 
culation of the single-particle spectral function using the 
Chebyshev polynomial method is fast, efficient and direct 
with less memory consuming, superior to direct diagonal- 
ization [generally 0(N 3 )] of the BdG matrix. 



III. NUMERICAL RESULTS 

In our work, the parameters are chosen as below: t = 1 
(as the unit of energy), t' = —0.3, fx = —0.83, and N = 
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FIG. 1: Upper-left panel: Superfluid density(SFD) and the long-range phase correlation function s(l4, 0) vs T for three different 
interactions. Black disks (SFD, V = 1.2), green square (SFD, V = 2.4), blue diamond (SFD, V = 4.0), orange down-triangle 
(S(14,0), V = 1.2), and red triangle (S(14,0), V = 4.0). T c is around 0.05, 0.09, 0.105 respectively. Upper-right panel: The 
d-wave vortex number n v vs T for three different interactions. For each interaction, there is a region that vortices gush abruptly. 
Black disks(U = 1.2), red square(U = 2.4), and blue diamond(U = 4.0). Panels of the Lower row: Snap shots of the spatial 
distribution of the magnitude and phase of the d-wave order parameter A<j(i) on a 29 x 29 lattice (note the periodic boundary 
condition) at T = 0.06 (lower left) and T = 0.12 (lower right). The phase on each lattice site is represented by a blue arrow, 
while the magnitude is represented by the size of the arrow as well as the gray scale density. Also shown are the topological 
excitations with the red Q denoting the vortex with winding number 1 and green (^) denoting the antivortex with winding 
number — 1. 



L x L = 28 x 28. According to the above parameters 
the average electron number is approximately 0.9 and 
accordingly the hole doping is 0.1. For each temperature, 
the first 10 3 MC sweeps are dropped to equilibrate the 
system. 10 3 configurations are used as samples to get 
statistical average. Each MC sweep includes 2N 2 local 
updates to reduce the configuration correlation. By these 
arguments, the statistical error reduces to an acceptable 
level. Then, we calculate the spectral function directly 
with the help of the Chebyshev polynomials^! with the 
truncation M = 2048. 

In Fig. [ija), we show the superfluid density (SFD) 
D s /ire 2 (see Appendix) as a function of temperature 
for different pairing interactions. The SFD decreases 
as the temperature increases and displays an apparent 



drop indicating a phase transition for the three inter- 
action strengths, although SFD has a long tail above 
the transition temperature due to the finite size effect. 
By taking the leading-edge midpoint as T c , we have 
T c w 0.05,0.09,0.105 for V = 1.2, 2.4, 4, respectively 
Compared with the previous worle^J, we find that 
the transition temperature increases with the interaction 
strength at least for V < 4. This increase of T c with V 
for small to intermediate value of V is further supported 
if we further examine the long-range phase correlation^ 

5(L/2,0) = ^Eili < e*** e'^W >, where <j>f de- 
notes the phase of the pairing field A(z, i+x). The results 
are also shown in Fig.jlja), indicating that both the SFD 
and S(L/2, 0) are measures of the phase stiffness of the 
condensate. 
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Conventionally, the phase fluctuation scenario relates 
the normal to SC phase transition in underdoped HTCS 
to the KT-type phase transition. Although people have 
made lots of efforts, it is still unclear how this occurs. 
Generally, fluctuations of the phase degrees of freedom 
of the superconducting order parameter can be described 
by the 2D XY model, the studies of which have re- 
vealed that the proliferation and unbinding of the vortex- 
antivortex pairs above Tkt destroys the quasi-long-range 
phase coherence. Here we explore this aspect by observ- 
ing the vortex- type excitations in the phase field (pd(i) 
of the d-wave order parameter A<j(i) = \Ad(i)\e lVd ^\ 
The definition of Aj(i) is A<i(i) = [A(i, i + x) + A(i, i — 
x) — A(i,i + y) — A(i,i — y)]/4. The (anti)vortices are 
plaquette-centered topological defects of the phase field. 
The winding number or vorticity of the (anti)vortex is 
defined as the anticlockwise sum of the phase differ- 
ence around each plaquette of the square lattice (di- 
vided by 2ir). For each plaquette labeled by its lower- 
left corner i, its four vertices anticlockwisely are i\ = 
i,%2 = i + x,is = i + x + y,i4 = i + y. The phase 
difference between two NN sites for instance i-i and i\ 
is 9 2 s{i) = <Pd(h) - <pd(h) = Imlog(A d (i 2 )A^(ii)). 
Therefore the winding number around the plaquette i 

is W(i) = [02,1 (0 + #3,2 («) + #4,3 W + 01,4(O]Ar- «>(*) = 1 

or —1 represents a vortex- or antivortex-type topologi- 
cal defect around the plaquette i. The total number of 
(anti) vortices, i.e. n Vl which quantifies the phase fluctua- 
tion relevant to the topological excitations, is calculated 
according to n v = J2i$w(i),i, i- e - the total number of 
plaquettes with w(i) = 1 (note that without external 
magnetic field the antivortex number always equals to 
the vortex number). The temperature dependence of n v 
is shown in Figs.JlJb). One can observe the abrupt jump 
of n v at approximately the same temperature obtained 
from Fig. [lja) , giving further indication of the KT-type 
phase transition. For illustration, Figs.JlJc) and (d) show 
snapshots of the pairing fields recorded at T = 0.06 and 
T = 0.12 for V — 2.4. For these two temperatures the 
coherence length of the d-wave order parameter is small 
and the vortices can be clearly identified as shown in the 
figures. At the temperature well below Tkt, the vortex 
density is dilute and all vortices are bound together as 
vortex-antivortex pairs as shown in Fig.JIjc). Above Tkt, 
vortices gush and the unbinding of vortex-antivortex pair 
is clearly illustrated as shown in Fig. |ljd) . This behavior 
conforms to the characteristics of KT phase transition. 
Together with the observation that the crossover regions 
of the SFD and phase correlation reduplicate that of the 
vortices, we can argue that the SC phase transition is of 
the KT type. 

The single-particle excitation spectra of electrons mov- 
ing in a fluctuating pairing field is highly nontrivial. Fig- 
ure [2] shows the evolution of A(k,uj — 0) as a function 
of temperature. The momentum k is in the first Bril- 
louin zone and the energy uj = on the Fermi level. 
We set the chemical potential fj, — —0.83 intentionally 
to have more discrete momenta on the underlying Fermi 




FIG. 2: Temperature dependence of the spectral function at 
the Fermi energy A(k., cj = 0) with k in the first Brillouin 
zone (V = 2.4): From (a) to (d) T = 0.03,0.09,0.15,0.3. At 
T = 0.03, four sharp spectal peaks at four nodes; With the 
increasing temperature, A(k, 0) decreases at the nodes while 
piles up at the other k's of the underlying Fermi surface; At 
T = 0.3, the spectral distribution seems like a bowl. The gap 
of the antinode survives at the high temperature (For 28 x 28 
lattice, T > 0.3; For smaller size lattice, we find it dose not 
close even at T > 0.5). 



surface and the resulting electron number is around 0.9. 
At T — 0.03, four sharp spectral peaks right at the four 
gap nodes are clearly resolved in Fig. [2ja), which indi- 
cates that at temperatures well below Tkt the pair fluc- 
tuations are rather weak and the Fermi surface are ac- 
tually point like as in pure d-wave superconductors. At 
T = 0.09, the height of the peaks falls while their profile 
extends towards the antinodal direction, i.e. the spectral 
weight of other k points along the underlying Fermi sur- 
face increases. From Figure [2j we find that this pile-up 
effect of spectral weight at the vicinity of the underlying 
Fermi surface increases with temper ature, which is con- 
sistent with the ARPES observation d 24 * 25 ! as well as the 
theoretical picture^. 

Next, we show the energy distribution of the spectral 
function in Fig. [3jb) and (c) to examine the spectral gaps 
opened at different k's on the underlying Fermi surface at 
two different temperatures. The selected three k's are the 
node (ff > ff ), the wave vector (ff , f ) near the node and 
the antinode (ff^, n) as shown in Fig. [3^ a). At T = 0.05 
below the KT transition, A(k, ui) for k at node displays a 
sharp peak located at zero energy. Away from the node, 
we find that the spectral gap opens at the selected mo- 
mentum k = (ff , f ) closest to the node and increases 
to its largest value at the antinode, which conforms to 
the characteristic of d-wave superconducting gap func- 
tion. At higher temperature T = 0.09, we find that the 
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T/t to/t T/t 

FIG. 3: (a), The Fermi surface shown as red line in the first quarter of the Brillouin zone. Three selected k points lying on 
the Fermi surface are shown as large blue dots. The spectral function A(k,u) as a function of cj for this three k points with 
(b) T = 0.05 and (c) T = 0.09. (d) The Fermi-arc length versus temperature, (e) Temperature dependence of A(k — ki4,cj). 
Black solid Iine(T = 0.05), red dashed line(T = 0.08), blue dotted line(T = 0.09), and green dot-dashed line(T = 0.1). (f), 
A(k, uj = 0) as a function of temperature for the node point and ks. For all panels V = 2.4. 



spectral peak at the node is lowered in consistent with 
the observation of Fig. [5Jb). Moreover the spectral gap 
at k = (ff , §) is closed, while a spectral peak is piled 
up at the zero energy, which signals the formation of the 
so-called Fermi arc that extends from node as far as to 
k = (ff , f )■ To quantify the length of the Fermi arc, we 
examine the loss of the spectral weighP^due to the open- 
ing of the spectral gap L = [1 — A(k|M i^!^ |M=A) ], 
where the spectral gap A is measured as half the peak 
to peak separation of the spectral function. In the ideal 
case L = 1 means the opening of a full spectral gap; while 
L = identifies the closing of the gap, or in other words 
the formation of the Fermi arc. Here in analyzing our nu- 
merical results, we set both L = 0.1 and L = 0.15 as the 
threshold for the arc formation. The results are plotted 
in Fig. J3jd) where the variation of the length of Fermi 
arcs as a function of temperature is shown. For both pa- 
rameters, the arc length exhibits apparent rise near T c , 
which is consistent with the APRES measurement!^. In 
addition, this jump locates around the same tempera- 
ture where SFD, phase correlation, and the vortex den- 
sity changes most remarkably, indicating the importance 
of pair fluctuation in the formation of the Fermi arc. We 



cut the underlying Fermi surface between the node and 
antinode in the first quarter of the first Brillouin zone into 
20 equally spaced parts, with ko denoting the node and 
k2o the antinode. We examine the ki4 point in Fig. [3^e) . 
It is clearly shown that there is a continuous increase of 
the spectral weight at the Fermi level from T = 0.05 to 
T = 0.1, while the spectral gap shrinks as the tempera- 
ture increases, which can be explained by the increasing 
broadening effect due to the thermal pairing fluctuation. 

Now we report the second shift of the zero energy spec- 
tral weight. According to Figs. |3jb) and (c), we notice 
that the zero energy spectral weight at the node/antinode 
is always decreaing/increasing with temperature. How- 
ever for the k points near the node, the zero energy 
spectral weight first increases with temperature and then 
decreases above a temperature whose value depends on 
k as shown in Fig. |3jf). We call this phenomenon the 
second shift, which can be understood according to the 
theory of Berg and AltmarW Because of the Doppler- 
shift effect, the zero-energy spectral weight of the node 
is transferred to its neighboring k points and is gradu- 
ally exhausted nearby T c ; the neighboring points received 
the zero energy spectral weight from the node, and also 
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due to the same effect, shift their zero-energy spectral 
weight to their neighboring points. The higher the tem- 
perature is, the less they get and the more they shift. For 
high enough temperature, both the node and the neigh- 
boring points have approximately equal amount of the 
zero-energy spectral weight as shown in figure [3](f ) , and 
the zero-energy spectral weights of these points saturate 
and begin to decrease. 



IV. CONCLUSIONS 

In conclusion, we have carried out the classical Monte 
Carlo simulation of the 2D attractive Hubbard model. 
We have presented a local-update procedure based on the 
Matsubara Green's Function using Nambu-Gor'kov for- 
malism, which speeds up the exploring of the configura- 
tion space. We found that thermal fluctuations do contri- 
bution to the pile up of low-energy spectral weight on the 
underlying Fermi surface and the evolution of fermi arcs 
with temperature. The abrupt jump of the arc length is 
a qualitative result caused by the continuous piling up 
of zero energy spectral weight, and the second shift sug- 
gests that E. Berg and E. Altman's idea works better 
than simply thermally broadening effect. Finally, tak- 
ing superfluid density as the SC criteria is effective when 
interaction V is small. 

The work was supported by the Natural Science Foun- 
dation of China No. 10674179. 



V. APPENDIX: SUPERFLUID DENSITY 

The superfluid density can be obtained applying the 
linear response theory to the homogenous superconduct- 



ing stated Here we will use Gor'kov Green's function to 
express the superfluid weight. We start from the follow- 
ing formula ^ 



ire* 



Nh 2 



(-K x ) - Il xx (q x = 0, q y -> 0,ifi m = 0), 



(37) 

where D s represents the superfluid weight and measures 
the ratio of superfluid density to mass. Here K x denotes 
the electron kinetic energy along x direction and 



K x 



(<L<**) _ oY - Gn(k^„) e -"° 4 



( 38 ) 

to£ = (d 2 E] si /dk x ) 1 the electron effective mass with 
£k = — 2£(cos k x + cos k y ) — 4i' cos k x cos k y the electron 
dispersion. Tl xx is the current-current correlation func- 
tion defined in momentum and imaginary-time space 



H xx (q,r) = ^<T r j x (q,T)i(-q,0)), (39) 



where j x (q) 



- i P W2 



E kCT w k+ q /2 c L c k+qa the current- 



density operator, and v%. — de-^/dk x denotes the group 
velocity of electron. Performing Fourier transform with 
respect to imaginary time, we have 



n ra (q,iO„)= / dTe^ T U xx {q, T ), 
Jo 



(40) 



Combining Eq. ( 37 ) , ( 38 1 and ( 40 1 followed by straight 



forward derivation, we have the equation for superfluid 
density expressed using the Gor'kov Green's function, 



rrr 



ire* 



Gn(k,iw„)e i 



Nh 2 ^ 

k.cj™ 



mi 



Nh 2 /3 ^ 



(< +q/2 ) tr[G(k , iui n )G(k + q, ioj n )] 



(41) 



r 



Then we will use the above formula to pair-fluctuating 
superconductors, whose superfluid weight is given by 



J DADAe-^ A A ) 



(42) 



where, D S (A, A) denotes the superfluid weight for a cer- 
tain configuration A. Considering that A is spatially 
inhomogeneous, we should first perform Fourier trans- 
form on the real-space Gor'kov Green's function which 



has been obtained and updated during the random walk 
through the configuration space, 



G(k, iu n ; A) = 1 £ G{i, j,iu n ; A)e ik < i ~^ 



(43) 



After the transformation, Eq. ( 43 1 is inserted into 



Eq. (41 1 to calculate the superfluid density corresponding 



to one configuration and then using Eq. ( 42 1 to obtain 
the statistical average over configurations. 
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